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ABSTRACT 

In this paper, we study small- iV gravitational dynamics involving up to six objects. We 
perform a large suite of numerical scattering experiments involving single, binary, and 
triple stars. This is done using the FEWBODY numerical scattering code, which we 
have upgraded to treat encounters involving triple stars. We focus on outcomes that 
result in direct physical collisions between stars, within the low angular momentum 
and high absolute orbital energy regime. The dependence of the collision probability 
on the number of objects involved in the interaction, N , is found for fixed total en- 
ergy and angular momentum. Our results are consistent with a collision probability 
that increases approximately as N'^ . Interestingly, this is also what is expected from 
the mean free path approximation in the limit of very large N . A more thorough 
exploration of parameter space will be required in future studies to fully explore this 
potentially intriguing connection. This study is meant as a first step in an on-going 
effort to extend our understanding of small- collisional dynamics beyond the three- 
and four-body problems and into the realm of larger-A. 

Key words: gravitation - (stars:) binaries (including multiples): close - methods: 
statistical - celestial mechanics - scattering. 



1 INTRODUCTION 

The gravitational three-body pro blem was firs t studied by 
Sir Isaac Newton in liis Principia (|Newtonlll686l ). After solv- 
ing the two-body problem, Newton boldly added a third 
body into the mix and attempted to create a similar math- 
ematical formalism to describe the motion of three celestial 
bodies under their mutual gravitational attraction. The so- 
lution evaded Newton until his end, leaving the problem 
at t he mercy o f his disciple s. It was later taken up by Eu- 



ler jEuleilll772h . Lagrange (lLagrangelll8l"l|). Jacobi (|jacobil 



118361), Poincare jPoincardI 18921 ). Hill (e.g. lHilllll878l ') , Henon 
(e.g. lHenonlll969i r and a host of ot hers, often with a focus 
on th e Earth-Moon-Sun system (e.g. IValtonen fc KarttunenI 
l2006l ). 

Despite the considerable efforts of numerous re- 
searchers, a simple analytic solution to the general t hree- 
body problem has never been found. ISundmanI (|l912') dis- 
covered a uniformly convergent infinite series involving fa- 
miliar functions that technically solves the three-body prob- 
lem. However, in order to attain a reasonable level of accu- 
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racy, t he solution must contain on the order of iq^.ooo.ooo 
terms IValtonen fc Karttuncn 2006). More practical solu- 
tions require a number of simplifying assump tions to mak e 
the general three-body problem tractable (e.g. lHenonlll969l ). 
As a result, the most useful analytic solutions tend to be 
solely applicable to a very narrow subset of the total allowed 
parameter space. 

The introduction of computers within the last few 
decades revolutionalized our understanding of the three- 
body problem. This allowed for the direct integration of 
the equations of motion for each particle over small time 
steps, incrementally moving each body forward in an itera- 
tive fashion. Despite the fact that this approach completely 
transformed our tool-set for studying the three-body prob- 
lem, it too has its short-comings. For example, the times re- 
quired to run the simulations to completion are often quite 
long. The precise definition of a "complete encounter" can 
often be ambiguous as well. Quasi-stable configurations can 
form that remain bound for millions or even billions of years 
before eventually breaking apart (e.g. iMikkola fc ValtonenI 
1986). There is also the issue of errors in computing the 
trajectories of the particles in position-space which are in- 
troduced at each time-step. These arise as a result of mov- 
ing one body forward at a time, instead of moving all bod- 
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ies simultaneously. Such errors can be minimized by taking 
suitably short time-steps. However, this comes at the often 
considerable cost of simulation run-time. 

When only three bodies are involved, qualitative de- 
scriptions of the outcomes of dynamical interactions are of- 



changes, fly-by s . and collis i ons 



Mikkolal 1983 



19841 : IHu^ 11984 : 



e. IHut & Bahcalll 


1982 


McMillan & Hut 


1996 



Fregeau et al.l 20041 ). However, the number of possible out- 



comes quickly increases with increasing TV, where A'^ is 
the total number of objects involved in the interaction 
l|Leigh &: Sills! 1201 ll ) . This complicates descriptions of the 
interactions, and introduces a considerable challenge in de- 
veloping a physical understanding of the evolution of the 
system. For example, nearly 100 generic outcomes are pos- 
sible for encounters involving six objects. Additional bod- 
ies not only increase the parameter space to be explored, 
they also increase integration times for computer simu- 
lations. Consequently, most previous numerical scattering 
studies considered only single-bina ry and, to a lesser ex- 
tent, binary-binary encounters (e.g. lHeggielll975l: IHu3 
McMillanI Il98^: ISigurdsson fc PhinnevI Il993l: iDavied 



1983 



1995 



Bacon. Sigurdsson fc Daviesll 19961 : Fregeau et al.l 2004 ). 

Many of the numerical scattering studies conducted to 
date considered equal-mass particles, and devoted their at- 
tention to studying the effects of va rying the initial rela - 
tive velocity or impact pa rameter (e.g. lHut fc Bahcallll 19831 ). 
For example, iHuj (|l983h found analytic approximations for 
high-velocity encounters, and provided simple formulae for 
the corresponding cross-sections for exchanges and ioniza- 
tions to occur. Similar analytic formulae were derived by 
I Mikkolal (Il984l ) for encounters involving two binaries having 
unequal orbital energies but equal masses. 

An extensive analys is that considered unequal mas s par- 
ticles was conducted bv ISigurdsson fc PhinnevI (|l993l ). who 
studied the effects of dynamics on the stellar and remnant 
populations in a globular cluster. A number of other scat- 
tering experiments were designed purely to study the forma- 
tion of different types of stella r exotica in glo bular clusters, 
including blue stragglers (e.g. 'Leonard!'l989l), low-mass x- 
ray binaries ( e.g. [S igurdsson fc Phinnev 1995), cataclysmic 
varia bles (e.g. [ ivanova et al.ll2006 ). and millisecond pulsars 
(e.g. Ilvanova et al. 20081 ). Manv of these also considered a 
range of particle masses. The effects of general relativity 
have also been studied in the c ontext of the three- and four - 
bod y problems. For example, iMikkola fc ValtonenI (|199G| ) 
and IValtonen et al.l l|l994h considered interactions involv- 
ing binary super-massive black holes during the mergers of 
galaxies, and identified several important trends arising from 
energy losses due to gravitational radiation. 

In this paper, we study gravitational interactions in- 
volving up to six objects. Our goal is to better understand 
how the outcome of an encounter depends on the number 
of interacting objects. This is a first step toward extending 
techniques developed for the three-body problem, where the 
vast majority of research efforts have thus far been concen- 
trated, to treat larger-A'^ encounters. To this end, we per- 
form > lO'^ numerical scattering experiments involving sin- 
gle, binary and triple star systems. The number of possible 
encounter outcomes is a steeply increasing function of A*'. 
This presents a considerable challenge when trying to draw 
parallels between encounters involving different numbers of 



objects. To minimize this issue, we focus only on outcomes 
resulting in direct physical collisions, which we consider to 
have occurred when the radii of any two stars overlap. 

In Section (2] we describe the set-up for our numerical 
scattering experiments, including the range of initial condi- 
tions considered. Additionally, we develop an analytic for- 
mula for the collision probability as a function of A'', that 
is based on a simple analytic model originally derived for 
1-1-2 encounters. In Section |31 we present the results of our 
analysis of this large suite of single-binary (1-1-2), binary- 
binary (2-f2), single-triple (1-1-3), binary-triple (2-1-3), and 
triple-triple (3-1-3) scattering experiments. Here, we fit the 
analytic formula to the results from these numerical scat- 
tering experiments, thereby deriving the A*'- dependence of 
the collision probability. The implications of our analysis 
for small- A'^ coUisional dynamics are discussed in Section [l] 
Concluding remarks are given in Section (5) Finally, in an 
appendix, we present a formalism for creating schematic di- 
agrams that quantitatively depict the evolution of an in- 
teraction in energy-space, and briefly discuss some of their 
possible applications. 



2 METHOD 

In this section, we describe the details of our numerical scat- 
tering experiments, and present the general form for the col- 
lision probability to which our results will be compared in 
Section |31 



2.1 Scattering Experiments 

We calculate the outcomes of a series of single-binary, 
binary-binary, single-triple, binary-triple, and triple-triple 
encounters using the FEWBOdM3 numerical scattering 
code. The code integrates the usual A'-body equations in 
configuration- (i.e. position-) space in order to advance the 
system forward in time. For details concerning the adap- 
tive integration, classification t echniques, etc. used b y FEW- 
BODY, we refer the reader to lFregeau et al.l (|2004l ). 

We adapted the FEWBODY code to handle all encoun- 
ters involving singles, binaries, and triples. (Specifically we 
created additional subroutines to simulate 1-1-3 and 3+3 en- 
counterfl; codes to simulate encounters between binaries 
and singles only, as well as a 2-1-3 encounter code, were pre- 
viously availabl e in the FEWBODY package.) We use the 
same criteria as lFregeau et al.l (|2004l ) to decide when a given 
encounter is complete, defined as the point at which the sep- 
arately bound hierarchies that make up the system are no 
longer interacting with each other or evolving internally. 

To perform physical collisions between stars, FEW- 
BODY uses the "sticky-star" approximation. This treats 
stars as rigid spheres with radii equal to their stellar radii. 
A physical collision is assumed to occur when the radii of 
the stars overlap. When this happens, the stars are merged 
assuming conservation of linear momentum and no mass 



^ for the source code, see http: / /fewbody.sourceforge.net" 
^ The authors are happy to provide these additional subroutines 
to users of FEWBODY upon request. They will be made avail- 
able on the FEWBODY homepage shortly, along with the latest 
version of the source code (FEWBODY 0.27). 
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loss. This does not account for tidal effects, which could sig;- 
nifica ntly increase the collision probability (e.g. iMcMillanI 
1 19861 ). but are beyond the scope of this work. For this rea- 
son, we consider the collision probabilities presented in this 
paper to be lower limits for the true values. 

Previous scattering experiments have shown that the 
total energy and angular momentum are the most impor- 
tant para meters in deciding the outco mes of l-f2 interac- 
tions (e.g. IValtonen fc Karttunenll200(j 'l. Therefore, we will 
fix the total energy and angular momentum when comparing 
between encounters involving different numbers of objects. 
By fixing these quantities, we aim to remove the depen- 
dences of the encounter outcomes on energy and angular 
momentum (when comparing between different encounter 
types), thereby normalizing the comparisons to reveal the 
A'^-dependence of the collision probability. We consider sev- 
eral different combinations of the total energy and angular 
momentum, which we henceforth refer to as Runs (see be- 
low). For each of these combinations or Runs, both the total 
energy and angular momentum are always chosen to be the 
same to within a factor of ~ 2 for every type of encounter. 

We focus on three primary Runs (labeled Runs 1, 2, and 
3 in Table [1} in Section [3] However, in order to check the 
robustness of our results, we also perform eight additional 
Runs with similar total energy and angular momentum as 
adopted in Runs 1, 2, and/or 3. We selected the following 
additional combinations of semi-major axes for the two or- 
bits used in each of these additional Runs: 0.5 AU and 3.5 
AU, 0.5 AU and 6.5 AU, 0.5 AU and 8 AU, 0.5 AU and 
10 AU, 1 AU and 8.5 AU, 1 AU and 11.5 AU, 1 AU and 
15 AU, 2 AU and 15 AU. For all three primary Runs, we 
perform a total of 800,000 numerical scattering experiments 
for each of the different encounter types, whereas this num- 
ber is reduced to 200,000 for the other Runs. In total, this 
amounts to 12 million simulations. This number is chosen to 
be a balance between statistical significance and simulation 
run-time, which can be quite long for simulations involving 
five or more objects. 

As we design our experiment such that the total angular 
momentum is roughly the same for every encounter type 
within a given Run, we choose the sum of the semi-major 
axes of the two objects involved in the interaction (which 
we will call the geometnc cross-section) to be equivalent 
to within a factor of ~ 2 for every encounter type. These 
cross-sections are determined by the initial semi-major axes 
of any binaries and/or triples involved in the interaction, 
since the physical radii of the stars are small in comparison. 
(For triples, the cross section is primarily determined by the 
semi-major axis of the outer orbit.) Specifically, we choose 
the semi-major axes shown in Table [T] A semi-colon is used 
to separate different objects, whereas a comma is used to 
separate the orbits within triples. Parantheses are used to 
enclose the semi-major axes of triples, with the smaller of 
the two semi-major axes always corresponding to the inner 
binary. 

We assume equal masses of 1 Mq and stellar radii of 
1 R0 for all stars, and circular orbits for all binaries and 
triples. This is done to make our exploration of the rele- 
vant parameter space more tractable. However we expect 
the collision probability to in general depend on the mass 
ratio s and orbital eccentricties (e.g. Sig urdsson fc PhinnevI 
1 19931 ). Unequal masses and non-zero eccentricities are be- 



yond the scope of the present paper. However we intend to 
address these issues in future work. The ratio between the 
inner and outer semi-major axes is chosen to be ^ 7 for 
all triples, since this roughly corresponds to the stability 
boundary for the internal evolution of triples composed of 
equal mass-particles with initially zero-eccentricity inner or- 
bits and moderate-eccentricity outer orbits (iMardling 200 1| ). 
The initial angle of inclination between the inner and outer 
orbits of all triples, in addition to their initial relative phases, 
are chosen at random. 

For each Run, we populate a grid of scattering experi- 
ments varying only the relative velocity at infinity and the 
impact parameter. Specifically, we select values for the rela- 
tive velocity at infinity from to 1.1 in equally-spaced inter- 
vals of 0.004, in units of the critical velocity Vcrit (defined as 
the relative velocity that gives a total energy of zero for the 
encounter). We select values for the impact parameter from 
to 7 in equally-spaced intervals of 1, in units of the geomet- 
ric cross-section for collision (e.g. the semi-major axis of the 
binary for a 1-1-2 encounter, the sum of the semi-major axis 
of the outer orbit of the triple and that of the binary in a 2-|-3 
encounter, etc.). In this way, for a given impact parameter, 
we ensure that the range in both the total angular momen- 
tum and the total energy are equal for every encounter type 
to within a factor of ~ 2. Finally, at each point in this grid, 
we perform multiple scattering experiments randomizing all 
angles in the encounter, including the angles at impact be- 
tween the orbital planes of any binaries and triples involved 
in the encounters. 



2.2 Collision Probability 

In this section, we present a general functional form for the 
collision probability. We begin by comparison to the 1-1-2 
encou nters studied in Section 8.6 of IValtonen fc KarttunerJ 
(200^, and then generalize this formula to A'^ > 3. The best- 
fitting values for all free parameters in this analytic collision 
probability formula are then found in Section [3] through fits 
to the results of our numerical scattering experiments. 



2.2.1 Collision Probability for N = 3 

The simple analytic model we use to guide our choice for 
the collision probability was originally found for resonant 
1-1-2 i nteractions. The model i s descr ibed in detail in Section 
8.6 of IValtonen fc KarttunenI l|2006l ). to which we refer the 
reader for its full derivation. It is based on the assumption 
that the encounter will evolve via a series of successive ejec- 
tions, eventually leading to the escape of one of the stars 
from the three-body system. We use the term ejection to 
describe a scenario in which one object ends up at a great 
distance from the others before falling back to resume the in- 
teraction. Some of these ejections are of considerably longer 
duration than others. Every ejection counts, however, since 
they all produce a temporary binary independent of the du- 
ration of the ejection. The temporary binary experiences 
close two-body encounters at the pericentre q of every orbit 
(e.g. ISzebehelvl Il967l : ISaslawl Il974l : IValtonen fc KarttunenI 
I2OO6I ). 

For a 1+2 interaction, the probability that two stars 
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Table 1. Initial semi-major axes of all binaries and triples for every Runs 1, 2, and 3 



Encounter Type Run 1 Run 2 Run 3 

(in AU) (in AU) (in AU) 

1+2 10.0 5.0 7.0 

2+2 1.0; 10.0 0.5; 5.0 1.0; 7.0 

1+3 (1.0, 10.0) (0.5, 5.0) (1.0, 7.0) 

2+3 10.0; (1.0, 10.0) 5.0; (0.5, 5.0) 7.0; (1.0, 7.0) 

3+3 (1.0, 10.0); (1.0, 10.0) (0.5, 5.0); (0.5, 5.0) (1.0, 7.0); (1.0, 7.0) 



will pass within a distance q from each other can be approx- 
imated by: 



Pcoiiiq) ~ 240C(L)g/ao, 



(1) 



where ao is the initial semi-major axis of the binary, and we 
require q <^ ao- The factor C{L) is needed to account for the 
fact that the lifetime of t he system depends on the total an- 
gular momentum L (e.g. IValtoneiilll974l : lAnosova" fc Orlovl 
Il986l ). If we take q to be the physical sizes of the objects 
involved in the encounter, then Equation [T] gives the colli- 
sion probability for a resonant 1+2 interaction. It has been 
shown to give good agr eement with the results of num erical 
scattering experiments (jValtonen fc Karttunen|[2006l ). 

2.2.2 Collision Probability for N > S 

We write the collision probability for an encounter involving 
A*' > 3 objects as: 

Pcoii{q)^aN''C{N,L)q/ao, 



where 



CiN,L) 



(2) 



(3) 



and a, /3, 5, and 7 are all constants. We take ao to be the 
semi-major axis of the shortest-period orbit involved in the 
interaction, and C{N,L) is a function of both the total an- 
gular momentum L and the number of objects N. As we will 
show in Section [3l Equation [2] gives excellent agreement to 
the results of our numerical scattering experiments. Below, 
we justify our choice for this functional form for the collision 
probability. 

First, we are interested in quantifying the dependence 
of the collision probability on the number of objects N in- 
volved in the interaction. This is accounted for in Equation[2] 
by a general power-law dependence on N (i.e. /?). Second, 
based on the results of previous numerical scattering exper- 
iments, we expect that an encounter outcome will depend 
both on the total energy and the total angular momentum. 
Therefore, we expect that these two quantities should ap- 
pear somewhere in the equation for the collision probability. 
These dependences are included in Equation[2]via the terms 
C{N,L) and q/ao. The first term C{N,L) is a function of 
the total angular momentum L, whereas the second term is 
roughly proportional to the total energy of the encounter. 
This last point follows from the fact that we have defined ao 
to be the semi-major axis of the shortest-period orbit, and 
it is this orbit that has the largest absolute orbital energy. 
Moreover, previous numerical scattering experiments of 1+2 
and 2+2 encounters have shown that the components of the 
shortest-period binary involved in the interaction are the 



most like ly to experien ce a collision during an encounter 
(e.g. Fre geau et aLl |2004). 

We allow for a possible A-dependence in the function 
C{N,L), since previous numerical scattering experiments 
performed to constrain this function considered only 1+2 
interactions. Therefore, we do not know whether the total 
number of stars involved in the interaction will affect its life- 
time, and play a role in determining the function C(A, L). 
As we will show in Section |3] the specific functional form we 
have adopted for C{N, L) in Equation [3] is needed to ensure 
that the correct agreement with the results of our numerical 
scattering experiments persists as we move to larger total 
angular momentum. 

It is important to note that Equation [2] does not ap- 
ply to 1+2 encounters. Therefore, we do not include them 
when finding the best-fit parameters. This is because the 
conditions imposed by our assumptions (e.g. equal masses 
for all stars, only circular orbits, etc.) are such that 1+2 
interactions cannot always be made to fit into our normal- 
ization for comparing between the different encounter types 
without significantly modifying the initial parameters of the 
encounter. Specifically, 1+2 encounters initially involve only 
a single bound orbit (via the binary) . However, all other en- 
counter types initially involve multiple orbits, which affords 
us additional free parameters. These can be adjusted to get 
the total energy and angular momentum to within our re- 
quired factor of 2, and therefore define a suitable normal- 
ization between encounter types. This is not always possible 
within the confines of our assumptions when 1+2 encoun- 
ters are also included. We will come back to this issue in 
Section [l] 



3 RESULTS 

The collision probability is calculated from the output of our 
numerical scattering experiments for a given encounter type 
and a given Run as: 

Ncou 



Pr.nll — 



N 



(4) 



where Ncoii is the total number of encounters that result 
in a direct physical collision, and N is the total number 
of encounters performed. The uncertainty for the collision 
probability is calculated using Poisson statistics according 
to: 



N 



(5) 



Technically, this uncertainty should include scattering ex- 
periments that result in unresolved outcomes (Hut 1983|). 
However, we find that the number of unresolved outcomes 
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is sufficiently small that it does not significantly contribute 
to APcoii, and we do not include it in its calculation. 

We show these collision probabilities as a function of 
the total angular momentum for 1+2, 2+2, 1+3, 2+3, and 
3+3 encounters in Figures [U [5] and [3] for Runs 1, 2, and 
3, respectively. The total angular momentum, shown on the 
X-axis, is provided in units of Mq cTgeoVcrit, where ag^o is 
the geometrical cross-section for collision in AU (which is 
determined by the semi-major axes of the orbits going into 
the interaction). 

We then fit Equation[2]to each of our individual Runs in 
order to derive the best-fitting values for the free parameters 
a, P, 5, and 7, that correctly predict the collision probability 
simultaneously for the 2+2, 1+3, 2+3, and 3+3 encounters. 
We use the Levenberg-Marquardt least-squares fitting tech- 
nique implemented in the MPFIT code (Mar kwardt 20081 : 
lMorelll978l ) to derive these best-fitting values. In order to ob- 
tain realistic uncertainties on these fit parameters, we scale 
the uncertainties on our measurements by a constant factor 
to make the reduced value equal to unity. 

The solid lines in Figures [Tl [2] and [3] show our fits to 
Equation [2] The best-fitting parameters are shown in Ta- 
ble [5] for all three Runs, along with their corresponding un- 
certainties. 

The mean power-law index for all 11 Runs is ;3 = 
2.40 ± 0.12. We see weak evidence for an increase in our /3 
values with increasing geometric cross-section, with a Pear- 
son's linear correlation coefficient of 0.73. We suspect that 
this weak trend is due to our normalization beginning to 
break down at large impact parameters, which can amplify 
small differences in the geometric cross-sections for different 
encounters types. The best-fitting power-law index on A'' is 
consistent with a value of /? = 2 to within roughly three 
standard deviations for all Runs. Moreover, a chi-squared 
test between all 11 of our j3 values and a constant power- 
law index of /3 = 2 returns no significant distinction, with 
a reduced chi-squared value of 1.73. Therefore, we conclude 
that our results are consistent with a collision probability 
that increases approximately as N"^ . 

Additionally, we find an exponential drop-off in the col- 
lision probability with increasing angular momentum that 
is steeper for larger N . Once again, this is the case for all 
Runs. We will return to these two features of the collision 
probability relation in Section H) 

We find some small disagreement between the collision 
probabilities for 2+2 and 1+3 encounters for Run 3. This can 
be understood as follows. For 1+3 encounters, the geometric 
cross-section for collision is equal to the semi-major axis of 
the outer orbit of the triple (since the radius of the single 
star is negligible in comparison). For 2+2 encounters, how- 
ever, the geometric cross-section is equal to the sum of the 
semi-major axes of the two binaries. Therefore, within the 
framework of our experimental set-up, the approximation 
that the geometric cross-sections for 1+3 and 2+2 encoun- 
ters will be equal only holds provided one of the binaries 
in the 2+2 case has a much greater orbital separation than 
the other binary. This assumption is a good approximation 
for Runs 1 and 2, however it begins to break down for Run 
3 with increasing impact parameter (i.e. increasing angular 
momentum in Figure [3]). 




100 200 300 

Total Angular Momentum 



Figure 1. The probability of a single collision occurring as a 
function of the total angular momentum is shown for Run 1 (see 
Tabic [TJ for every type of encounter. The open stars correspond 
to 3+3 encounters, the solid triangles to 2+3 encounters, the open 
circles to 1+3 encounters, the solid circles to 2+2 encounters, and 
the crosses to 1+2 encounters. The solid lines show our best-fits 
to the data using Equation [2] The top line corresponds to the 
case N = 6, the line below to TV = 5, and the bottom line to 
A'^ = 4. We show the scaled uncertainties here (and in Figures |2] 
and|3]l, as discussed in Section [3] 



4 DISCUSSION 

Within the angular momentum and energy regime stud- 
ied here, we find that the probability of a direct physical 
collision occurring during an interaction increases roughly 
as . One interpretation of this result can be under- 
stood as follows. Previous numerical scattering experiments 
of 1+2 interactions have shown that the collision proba- 
bility is proportional to the average number of close ap- 
proaches experienced by the system per crossing time, mul- 
tiplied by the number of cro ssing times survived through 
jValtonen fc KarttunenI I2OO6I) . Therefore, this depen- 
dence may arise physically from an dependence in the 
total number of close approaches experienced by the system. 
This can occur in one of (at least) two ways. Either the num- 
ber of close approaches per crossing time scales as while 
the number of crossing times survived through is constant, 
or the number of crossing times the system survives through 
(until the time when the first collision occurs) scales as 
while the number of close approaches per crossing time is 
constant. 

We find that the average number of crossing times at 
the time of the first collision is the same to within roughly 
a factor of ~ 2 for all encounter types, and that there is no 
trend with N. Therefore, it is unlikely that the latter sce- 
nario described above is responsible for producing the 
dependence. This suggests that perhaps the A^ dependence 



© 2011 RAS, MNRAS 000,[lHIi] 



6 Leigh & Geller 



Table 2. Best-Fit Parameters for Equation [2] 



Parameter 


Run 1 


Run 2 


Run 3 


a 


1.10 ± 0.20 


1.06 ± 0.19 


1.73 ± 0.32 




2.39 ± 0.11 


2.15 ± 0.11 


2.25 ± 0.11 


5 


-6.47e-2 ± 4.40e-3 


-4.60e-2 ± 3.21e-3 


-1.15e-l ± 9.97e-3 


7 


1.33e-l ± 1.56e-2 


1.33e-l ± 1.56e-2 


2.31e-l ± 2.33e-2 




200 400 

Total Angular Momentum 




50 100 150 200 

Total Angular Momentum 



Figure 2. The probability of a single collision occurring as a 
function of the total angular momentum is shown for Run 2 (see 
Table [TJ for every type of encounter. The symbols representing 
each encounter type, as well as the lines showing our best-fits to 
the data using Equation [2] for each encounter type, are the same 
as in Figure [T] 



Figure 3. The probability of a single collision occurring as a 
function of the total angular momentum is shown for Run 3 (see 
Table [TJ for every type of encounter. The symbols representing 
each encounter type, as well as the lines showing our best-fits to 
the data using Equation [2] for each encounter type, are the same 
as in Figure [T] 



in the collision probability may arise from a similar depen- 
dence in the number of close approaches per crossing time. 

However, in practice, this hypothesis is much more dif- 
ficult to test quantitatively. In particular it is not immedi- 
ately clear how to define a "close approach". We attempt 
to quantify this effect by examining animations of a handful 
of encounters in position-space, and counting the number of 
close approaches per crossing time by eye. It is clear that the 
number of close approaches per crossing time increases with 
increasing N . However, this method is far from adequate to 
determine the precise A^-dependence of this relation at any 
statistically significant level. 

In Appendix [X] we present a method that will im- 
prove upon this component of our analysis in future stud- 
ies. Specifically, we present a prescription for generating 
schematic diagrams that depict the evolution of an interac- 
tion in energy-space. As is explicitly shown in Appendix 1X1 
the advantage of this technique is to provide a straight- 
forward means of defining the criterion of "close approach" 
in terms of the fraction of the total energy exchanged be- 
tween two individual stars. This directly relates the defini- 
tion of close approach to the total energy and therefore the 



initial conditions of the encounter, which is not the case if 
this criterion is defined purely in terms of a distance. 

Also note that, in a sense, a collision could be inter- 
preted as a very strict definition of a "close approach" . Here 
multiple crossing times are required before this definition 
of close approach is satisfied. As we find that the number 
of crossing times until the time of first collision is roughly 
equivalent (to within a factor of « 2) for all encounter types, 
our result can potentially also be expressed as the number 
of "close approaches", defined in this manner, per crossing 
time scaling as A'^^ . We will explore a range of other criteria 
to define a close approach in a future study and investigate 
explicitly the dependence of the number of close approaches 
per crossing time on the total number of stars involved in 
the interaction. 

Intriguingly, the collision rate, as derived from the mean 
free path approximation (e.g. iLeonar also scales as 

A'"^. In the limit of very large A'^, we would expect our re- 
lation for the collision probability to agree with what is 
predicted from the mean free path approximation. On the 
other hand, in the limit of small- A'', it is not clear that the 
standard assumptions of the mean free path approximation 
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should still hold. For the angular momentum and energy 
regime studied here, the collision probability appears con- 
sistent with that of the mean free path approximation, at 
least for N = 4,5, 6. Further study is required to determine, 
e.g., how different energy and angular momentum regimes 
affect this relationship, what the minimum A'" is for which 
the standard mean free path approximations are still valid, 
etc. 

The second most noticable similarity between all Runs 
is an exponential decline in the collision probability with 
increasing angular momentum that is steeper for larger A*'. 
This effect is small, however, relative to the total drop in the 
collision probability as A*' decreases within a given Run (see 
Figures[T][2l and[3]). One possible explanation that is consis- 
tent with our results is that the number of close approaches 
per crossing time may decrease with increasing angular mo- 
mentum more steeply for larger encounters. 

Our results are applicable to a regime of low angular 
momentum and high absolute orbital energy. As the inte- 
gration time increases with increasing angular momentum, 
we focussed our attention on the low-angular momentum 
regime in this paper to maximize the number of simulations 
performed for each Run, and thereby increase the statistical 
significance of our results. 

We cannot probe lower angular momentum encounters 
without violating our assumptions. There are two reasons for 
this. First, the lower boundary for the long-term stability of 
triple systems corresponds to a ratio between t he inner and 
outer orbital separations > 7 (|Mardlingj l200ll ). Therefore, 
we cannot lower the total angular momentum by decreasing 
the semi-major axis of the outer orbit of the triple without 
simultaneously reducing the semi-major axis of its inner or- 
bit. This brings us to our second requirement, namely that 
q «C ao (where ao is the semi-major axis of the shortest- 
period orbit initially involved in the interaction and q is the 
physical sizes of the objects involved in the encounter). We 
performed two additional Runs with ao — 0.1 AU to test 
the limit of the assumption g ^ ao. In both cases, the re- 
sulting chi-squared values found using the best-fits for our 
free parameters in Equation [5] were significantly larger than 
we found for our other Runs. We interpret this as being due 
to a breakdown of the requirement g <C ao . Our results sug- 
gest that the assumption g ^ ao holds provided ao > lOOg. 
This is only a rough estimate for the lower limit for the ratio 
g/ao, and more work will be needed to better constrain its 
precise value. 

In future work, we intend to address the A^-dependence 
of the collision probability for higher angular momentum en- 
counters, as well as different mass-ratios and eccentricities. 
A non-circular orbit provides an additional free parameter 
that can be changed to affect the total angular momentum, 
but not the total energy. Therefore, the addition of a non- 
zero eccentricty should in principle allow us to include 1+2 
encounters in our estimates for the collision probability us- 
ing the same or similar normalization method (i.e. fixing 
the total angular momentum and energy) as adopted in this 
paper to compare between the different encounter types. 



5 SUMMARY 

In this paper, we perform a large suite of numerical scatter- 
ing experiments to study the probability of a direct phys- 
ical collision occurring during single-binary, binary-binary, 
single-triple, binary-triple, and triple-triple interactions. We 
quantify the dependence of the collision probability on the 
number of objects involved in the interaction A'^ for fixed 
total energy and angular momentum. 

Our results suggest that the collision probability in- 
creases approximately as A'^^, for A'^ — 4,5,6. This result 
is consistent with the hypothesis that the average number 
of close approaches per crossing time also scales as A''^, and 
we will fully investigate this possibility in a future study. 
Interestingly, this same A'^-dependence is predicted by the 
mean free path approximation in the limit of very large A'^. 
This similarity is rather intriguing, and further work inves- 
tigating the connection between these two regimes in A'^ will 
be highly valuable for our understanding of collisional dy- 
namics in the realm of not-so-small- A'^. 



APPENDIX A: SCHEMATIC 
REPRESENTATION OF LIOUVILLE'S 
THEOREM 

In this appendix, we describe the construction of schematic 
diagrams that quantitatively describe the evolution of an 
interaction in energy-space. 

Al Motivation 

Previous numerical scattering experiments for 1-1-2 interac- 
tions have demonstrated that the probability for a dynami- 
cal encounter to result in a direct collision is approximately 
proportional to the number of close approaches per crossing 
time multiplied by th e number of crossing times th e system 
survives through (e.g. IValtonen fc Karttunenll200d ). This is 
the basis for the analytic model described in Section 12.2.11 
to predict the collision probability during resonant 1-1-2 en- 
counters, and we have extended this formalism to also de- 
scribe encounters with A'^ > 3. We briefly discuss the difficul- 
ties in determining the precise number of close approaches 
during an encounter from an analysis of position-space in 
Section ID 

Here we present an improved method for determining 
the number of close approaches per crossing time by analyz- 
ing encounters in energy-space. There are three clear ben- 
efits to this technique, which we will expand upon below. 
First, defining the condition for a close approach in terms of 
energy rather than position maps each close approach to a 
nearly discrete change in energy, which is easily calculated. 
Second, the energy-space diagrams presented in this paper 
contain quantitative information pertaining to the relative 
energies of the objects during each stage of the encounter, 
whereas position-space diagrams do not (and are often very 
hard to interpret visually). Third, this method will poten- 
tially allow for a more general formalism to describe the 
collision probability during a dynamical encounter. 

Below, we outline our method and provide a few ex- 
ample scenarios. We conclude by describing the future work 
that will be enabled by this energy-based approach. 



© 2011 RAS, MNRAS 000,[lHIi] 



8 Leigh & Geller 



A2 Liouville's Theorem 

Liouville's Theorem is a key postulate in classical statisti- 
cal mechanics. It states that the time evolution of the dis- 
tribution function corresponding to an Hamiltoni an system 
is co nstant along any trajectory in phase space (|Liouvillel 
Il838t ). Said another way, the density of points represent- 
ing particles in x,p phase space is conserved as the system 
evolves in time, where x and p denote the 3-D position and 
momentum vectors, respectively, of the particles. This is a 
remarkable and powerful result that can be applied to any 
dynamically-evolving system provided the forces acting on 
the particles are conservative and differentiable. This last 
condition excludes collisions due to the sudden introduction 
of additional forces, such as tides, shocks, etc. 

As a system evolves dynamically, energy and angular 
momentum are exchanged between particles. This diffusion 
in energy- and angular momentum-space governs the trajec- 
tories of the particles through position- and velocity-space. 
Therefore, it is often the case that consideration of the evo- 
lution of a system in energy- and angular momentum-space, 
as opposed to position- and velocity-space, will more directly 
reveal the underlying physical processes that determine its 
outcome. As we will show in the subsequent sections, it 
follows that patterns in the dynamical evolution of an in- 
teraction are often most apparent in energy- and angular 
momentum-space. 

A3 Diagrams for N = 4 

We will begin by describing the rules for the construction 
of our schematic diagrams for the case A*' = 4. The total 
energy for a four-body encounter can be written: 

4 4 

where rrii is the mass of object i, Vi is the speed of object 
i with respect to the centre of mass of the system, and r^j 
is the distance separating objects i and j. We note that 
rij = lli — ijl, where r^ is the vector separating object i from 
the centre of mass of the system. We re- write Equation lAll 
in the form: 

i=l J = l;j7^i 

or 

4 

j=i 

Now, we can form a polygon by setting each of the an- 
gles equal to: 

a. = -(^) x 360°. (A4) 

(Note that the angle 360° here is a result of the more gen- 
eral formula 180° x (iV - 2), when = 4.) As the system 
evolves, the total energy E remains conserved but the in- 
dividual energy terms Ei will change. This is manifested 
visually via the transformation of our schematic diagrams 
with successive time-steps. The rules of our diagrams are 
such that the individual angles of the polygon change as the 
system evolves, while their sum remains constant. 



Bound objects are connected by solid lines, whereas ob- 
jects that are unbound are connected by dashed lines. Note 
that, if an object is unbound, then Equation IA4I suggests 
that the angle corresponding to its vertex will be negative. 
We will come back to this below. 

To illustrate our methodology, we show the evolution 
in energy-space of a 2+2 encounter in Figure lAll The cor- 
responding position-space diagram is shown in Figure IA2I 
The latter figure shows the evolution of a 2-1-2 encounter in 
position-space projected onto the x-y plane. In this interac- 
tion, two identical binaries that are composed of equal-mass 
(1 Mq) stars with semi- major axes of 5 AU approach from 
the left and right (we denote the initial time by t = to) to 
meet at the origin. In Figure IXTl Stars 1 and 2 comprise one 
of the binaries, and Stars 3 and 4 comprise the other. Both 
of these bound pairs are connected with solid lines, respec- 
tively, while all remaining pairs of stars are initially unbound 
from each other, and therefore connected via dashed lines. 
Shortly after the binaries meet at the origin in Figure IA21 
one star (Star 2) is ejected from the system, and is expelled 
toward the lower left of the diagram (t — ti). Star 2 is now 
an unbound single object, and hence is connected to all other 
stars via dashed lines in the top left inset of Figure [XT] The 
angle corresponding to its vertex is also now negative, since 
Star 2 has a positive total energy. 

The three remaining stars then undergo a resonant in- 
teraction as they drift toward the upper right of the diagram 
in Figure [X2I due to conservation of momentum {t = t2). At 
roughly {x,y) = (36,61), another star (Star 4) is ejected 
toward the lower right. The left-over binary (composed of 
Stars 1 and 3) leaves the diagram at the upper right, with 
a semi-major axis (and hence orbital energy) that is smaller 
than those of the two initial binaries (t — t^). The final state 
of the system is depicted in energy-space in the lower right 
inset of Figure lAll Only Stars 1 and 3 are connected via 
a solid line, since they are the only two stars that remain 
bound. The angles corresponding to the vertices for Stars 2 
and 4 are now negative, and the angles corresponding to the 
vertices for Stars 1 and 3 are both positive and obtuse (i.e. 
> 180°). Note that, although described in detail here, the 
specific evolution of the encounter is quite difficult to follow 
visually in position-space without the help of the previous 
text. 

As a further step, we depict this encounter in Figure [Xsl 
using th e formalism for maki ng scattering diagrams pre- 
sented in lHut fc Bahcall (|l983l ). In this plot, time increases 
from left to right, whereas vertical transitions denote the for- 
mation and/or destruction of temporary hierarchies in which 
two or more stars are in close proximity (i.e. strongly bound) 
to each other relative to all other stars. All four states of the 
system depicted in Figure [XT] are shown in Figure Esl The 
initial (i.e. t — to) and final (i.e. t = ts) states correspond 
to those shown at the far left and far right, respectively, 
of Figure IA3I The state shown at t = ti corresponds to a 
time shortly after Star 2 has been ejected from the system. 
Note that, within the triangle formed from the vertices of 
Stars 1, 3, and 4 at t = ti, all of the angles are comparable. 
This means that a temporary hierarchy has formed where 
Stars 1, 3, and 4 all have similar energies. This is represented 
in Figure IA3I where the first convergence of all three lines 
(for Stars 1, 3, and 4) occurs simultaneously. This scenario 
occurs once more (i.e. at t = 12) before the interaction is 
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t = t, 




t = t. 



t = t„ 



t = 



Figure Al. Schematic diagram showing the evolution in energy- 
space of a 2+2 interaction between identical binaries. The panels 
have been arranged in chronological order, starting at the lower 
left inset and rotating clockwise (i.e. to < ti < t2 < ^3) ■ A 
detailed description of this interaction has been provided in the 
text. The state of the system depicted in the inset corresponding 
to t = 42 shows one example of the many transformations of our 
diagram that occur between ejection events. This is due to the 
occurrence of several close approaches. 



eventually terminated by the escape of Star 4 from the sys- 
tem, leaving Stars 1 and 3 bound in a binary (which has 
a smaller orbital separation than the initial binaries going 
into the encounter). 

To summarize, after Star 2 is ejected, a three-body sys- 
tem remains that evolves via the formation of a temporary 
binary that mediates a series of successive ejections of the 
remaining single star (which can also be exchanged into and 
out of the temporary binary) . Eventually, the single star re- 
ceives a positive total energy and escapes from the system. 



A4 Future Work 

One primary application for this method is to determine 
the number of close approaches per crossing time for an 
encounter involving A*' stars. We outline this method here. 
A "close approach" can be defined using Equation IA4I to 
convert a minimum distance (i.e. q in Equation [T| into a 
minimum ang le0 If we require that the distance of "close 
approach" is comparable to the physical sizes of the objects 
(which are assumed to be small relative to the characteristic 
size of the interaction region and the initial semi-major axes 
of any orbits going into the encounter), then these events 



^ Note that in the point-particle limit a "collision" is considered 
to occur when two objects are within some minimum distance of 
each other. Therefore, our energy-space diagrams can be used to 
depict collisions. 
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Figure A2. Evolution of an 2+2 interaction in position-space 
projected onto the x-y plane. The binaries are composed of 1 
Mq stars, and have semi-major axes of 5 AU and circular orbits. 
The encounter involves an essentially parabolic collision between 
two identical binaries composed of equal mass stars, initially ap- 
proaching from the left and right. Distance is measured in units of 
the initial semi-major axes a of the binaries. When the encounter 
is over, two single stars and a binary emerge. 



correspond to times of near instantaneous and significant 
energy-exchange between two individual stars. This is be- 
cause the stars tend to be strongly accelerated or deceler- 
ated during very close approaches, so that significant energy 
and/or angular momentum can be exchanged. It follows that 
a "close approach" will appear as a near discrete transfor- 
mation of our energy-space diagrams, whereas the evolution 
appears more continuous and ch aotic in position - space . This 
is demonstrated in Figure 3 of iHut fc Bahcalj (| 19831 ) and 
also in Figure IA3I Therefore, individual close approaches 
will correspond to high values for dai/dt, where ai is the 
angle corresponding to star i in an energy-space diagram. 
Consequently, one can define a "close approach" as occur- 
ring when dai/dt > 77, where is a parameter that can be 
varied to specify the precise definition of a "close approach" . 
We intend to study the time-evolution of the quantity doi /dt 
in future work. 

As is evident from Figure IA2I the time evolution of 
the system is very difficult to observe in position-space, and 
counting the number of close approaches is essentially im- 
possible from this type of a diagram. In contrast, we can eas- 
ily count ten (to first order) close approaches in Figure IX3l 
To at least first order, each time a close approach occurs, 
one of the the lines in a scattering diagram makes a vertical 
transition. 
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Figure A3. Scattering diagram for the 2+2 interaction shown 
in position-space in Figure IX2l and in energy-space in Figure IXT] 
The time evolution of the interaction proceeds from left to right. 



This energy-space formalism is easily extended to depict 
the evolution of encounters involving more than 4 objects. 
(In principle, there is no limit to the number of objects that 
can be described with this method.) Furthermore, we can, in 
principal, use this formalism to re-write the functional form 
for the collision probability shown in Equation [2] in terms 
of only the number of objects A*", the total angular momen- 
tum, and energy. This is because our method re-defines the 
condition for a "close approach" in terms of energy, instead 
of distance q. The term g/ao in Equation [2] can thus be re- 
placed with a new term that depends directly on energy. 
This will further generalize our approach, and should be an 
improvement over the normalization used in this paper to 
compare between different encounter types within a given 
Run. 

In addition, the same method we have presented to de- 
pict the evolution of the system in energy-space can be modi- 
fied to depict its evolution in angular momentum-space. The 
prescription for this is analogous to the formalism presented 
in the previous sections for energy-space. In general, our 
technique can be used to depict the time evolution of any 
conserved quantity. 

We intend to use our diagrams in future studies to sim- 
ulate the dynamical evolution of small-A gravitational in- 
teractions in both energy- and angular momentum-space. 
This will be done in (accelerated) real time, so that we may 
simultaneously view the evolution of a system in position-, 
energy-, and angular momentum-space. This new tool will 
be a useful addition to studies with the over-arching goal of 
identifying trends in the evolution of dynamical interactions. 
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